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We  consider  a  simple  model  of  the  flux  in  a  rf  superconducting  quantum-interference  devic* 
(SQUID)  ring  subjected  to  an  external  periodic  magnetic  fleld.  The  dynamic  equation  describing 
the  flux  response  of  the  SQUID  is  solved  analytically  in  the  absence  of  damping  and  external  driving 
terms.  We  then  introduce  these  terms  as  small  perturbations,  and  construct,  for  this  system,  the 
Melnikov  function,  the  zeros  of  which  indicate  the  onset  of  homoclinic  behavior.  For  the  parame¬ 
ter  values  under  consideration,  excellent  agreement  is  obtained  between  our  theoretical  predictions 
and  numerical  calculations  of  the  stable  and  unstable  (i.e.,  time-reversed)  solution  manifolds.  A 
chaotic  attractor  is  shown  to  appear  somewhat  above  the  homoclinic  threshold. 


I.  INTRODUCTION 

In  this  paper  we  wish  to  consider  the  appearance  of 
homoclinic  instabilities  and  chaos  in  the  driven  rf  SQUID 
(superconducting  quantum-interference  device).  Since 
Poincare,'  it  has  been  known  that  under  perturbation,  the 
stable  and  unstable  manifolds  emanating  from  a  hyper¬ 
bolic  fixed  point  are  no  longer  identical  and  may  cross, 
giving  rise  to  an  infinite  number  of  intersections  (homo¬ 
clinic  points),  the  resulting  motion  being  so  complicated 
that  it  may  be  characterized  as  chaotic  (or  statistical). 
An  existence  theorem  of  Smale  and  Moser^  states  that  the 
motion  in  a  region  near  a  homoclinic  point  is 
homeomorphic  to  a  Markov  shift  map.  Thus,  in  this  re¬ 
gion,  the  test  of  the  wild  instability  of  the  motion  is  the 
presence  of  homoclinic  crossing.  Since  the  separatiix 
solutions  are  so  sensitive  to  perturbation,  a  simple 
theoretical  test  function  due  to  Melnikov^"*  may  be  used 
to  determine  the  presence  of  the  homoclinic  instability. 

This  simple  theoretical  analysis  has  been  applied  to  a 
number  of  driven  oscillators'*’’  and,  in  particular  to  the 
rf-driven  Josephson  junction,’""  which  constitutes  the 
primary  element  of  the  rf  SQUID  under  consideration  in 
this  work.  The  mechanical  analog  to  the  Josephson  junc¬ 
tion  is  a  driven  damped  pendulum  which  is  well  known, 
numerically,  to  exhibit  chaos  and  a  “strange”  attracting 
set.  Both  the  unperturbed  Josephson  Junction  and  the  rf 
SQUID  are  multistable  systems.  In  particular,  the  rf 
SQUID  exhibits  hysteresis  above  a  certain  threshold 
value  of  a  characteristic  parameter  [the  response  of  the 
device  below  its  hysteretic  threshold  to  a  general  pertur¬ 
bation  of  dc,  random,  and  periodic  components  has  been 
studied  by  one  of  the  authors”  (A.R.B.)].  In  this  work 
we  shall  concern  ourselves  with  the  operation  of  the  rf 
SQUID  in  the  hysteretic  regime.  Far  fewer  studies  con¬ 
cerning  the  appearance  of  chaos  have  been  done  on  the  rf 
SQUID  than  in  the  Josephson  junction  (for  a  review  of 
the  latter,  see  Ref.  8).  Dmitrenko  et  al.^^  report  an  exper¬ 


imental  broadband  amplitude  spectrum  in  the  hysteretic 
regime.  Smith  et  a/.'^  report  on  the  results  of  numerical 
simulations  suggesting  period  doubling  to  chaos  and 
qualitative  agreement  with  experiment.  Finally,  Ritala 
and  Salomaa*’  have  observed  the  transition  from  quasi¬ 
periodicity  to  chaos  via  the  Feigenbaum  period-doubling 
scenario  for  the  case  of  the  sinusoidally  driven  rf  SQUID. 
Their  work  represents  a  major  step  forward  in  our  under¬ 
standing  of  the  nonlinear  dynamics  associated  with  this 
system.  It  is  our  object  in  this  work  to  theoretically  pre¬ 
dict  the  onset  of  homoclinic  instability  by  means  of  the 
Melnikov  test  function,  and  then  to  compare  this  result 
with  a  numerical  calculation  of  the  manifold  crossing. 
This  latter  has  not  been  done  for  even  the  rf-driven 
Josephson  junction.  Comments  are  also  made  concerning 
the  appearance  of  a  strange  attractor  in  the  rf  SQUID 
above  its  homoclinic  threshold. 

In  its  simplest  form,  the  rf  SQUID  consists  of  a  single 
Josephson  junction  shorted  by  a  superconducting  loop 
having  an  inductance  L,  An  external  magnetic  field  pro¬ 
duces  a  geometrical  flux  in  the  loop  together  with  a 
circulating  supercurrent  /(t)=  — /j  sin(27r<I>/d>o),  where 
‘I>=<l>j-Fl,i  is  the  actual  flux  sensed  by  the  loop  in  the 
steady  state,  d>o  being  the  flux  quantum  (<l>o=/i/2e 
=2.07X10"*’  Wb).  The  flux  d>  in  the  SQUID  ring 
obeys  the  dynamical  equation'^ 


x  Pe 

— -i-T,x-l-x-f  — sin(2irx)=x^  , 


(1.1) 


where  the  sinusoidal  contribution  arises  from  the  Joseph¬ 
son  screening  current.  Here,  the  dot  denotes  the  time 
derivative,  X  =  <I>  /%,  =  /^q,  (o\=\/LC,TimL/R, 
and  P^^2,tLIi/%.  C  and  R  are  the  capacitance  and 
normal-state  resistance  of  the  loop,  /,  being  the  junction 
critical  current.  It  is  worth  pointing  out  that  the  quanti¬ 
ty  (l3g/2tr)(OQ  is  simply  the  plasma  frequency  oj]  of  the 
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junction;  Wy  is  the  frequency  of  low-amplitude  solutions 
of  the  Josephson-junction  equations  in  the  absence  of  an 
external  magnetic  flux.  The  parameter  P=Pf  /lit  deter¬ 
mines  the  hysteric  threshold  of  the  SQUID;  above  t:  criti¬ 
cal  value  p=Pc>  the  solutions  of  x  of  (1.1)  are  mul¬ 
tivalued.  For  the  undriven  (autonomous)  rf  SQUID  with 
negligible  damping,  one  may  case  (1.1)  in  the  form 


where  the  potential  U  (z)  is  given  by 

2 

+2irP(i)ls\n2  .  (1.3) 

Here,  for  later  convenience,  we  have  made  the  phase 
transformation  z=27rx —(ir/2),  x,  =Xq  being  an  external 
(background)  dc  signal.  One  may  readily  verify  via  sim¬ 
ple  graphical  analyses  that,  for  Xg =0,1,2, . . . ,  one  ob¬ 
tains  multivalued  solutions  to  (1.1)  above  a  critical  value 
)3j»0.7325.  For  the  cases  Xo=7,|,|-. ..  one  obtains 
Pe={2ir)~K  For  any  other  value  of  Xg,  the  critical  non¬ 
linearity  parameter  hes  between  the  two  limiting 
values  given  above.  Throughout  this  work  we  shall  set 
the  dc  driving  term  Xg  equal  to  zero;  the  potential  (1.3)  is 
then  symmetric  about  z  =  —-rr/2.  This  potential  has  been 
plotted  in  Fig.  1  for  p—2  and  6)g  =  1.  In  the  later  sections 
we  will  include  a  damping  force 

F^  =  -km,  (1.4) 

and  a  periodic  driving  force 

x,=F^  =  ^  sinM/-/g)]  ,  (1.5) 

in  (1.2). 

It  is  apparent  from  the  potential.  Fig.  1,  that  z,  and  Zj 
are  unstable  hyperbolic  fixed  points  and  that  five  elliptic 
fixed  points  occur  (for  the  range  of  z  values  used  in  this 
figure  in  general,  one  expects  an  infinite  number  of  fixed 
points).  In  this  work  we  will  focus  on  the  separatrix  ener¬ 


U(z)= 


al 


z-t-Y-27rxo 


z 


FIG.  1.  Potential  I/(2)vszfor(/?,«o,Xo)s(2, 1,0). 


gy  U[Zi)=U{z2)=U{z^)=U{z^)  and  thus  be  primarily 
concerned  with  the  region  of  the  two  lowest-lying  hyper¬ 
bolic  fixed  points.  The  steady  solutions  to  Eq.  (1.2)  are 
given  by 


— 27r/?cosZj  . 


(1.6) 


Simple  graphical  analysis  shows  that  for)3=s0.7325  a  new 
fixed  point  (other  than  z^  =  —Tr/2)  appears  at  Zj =7r.  For 
P  >0.1325  multistability  occurs  (as  pointed  out  earlier). 
The  character  of  these  new  fixed  points  may  be  analyzed 
by  linear  stability  analysis  with  the  eigenvalues  A.  of  the 
linearized  equations  determined  in  the  familiar  fashion.^ 
We  find 


A.=  ±£ag(l— 2irj8sinzj)'''^ ,  (1.7) 

from  which  it  may  be  seen,  for  P>Pc  ^0.1325,  that  the 
new  fixed  points  appear  (symmetric  about  z^  =  —it/2)  in 
“pairs”  of  hyperbolic  and  elliptic  points,  as  indicated  by 
the  potential  in  Fig.  1.  It  may  be  shown  that  the  hyper¬ 
bolic  point  Z|  in  Fig.  1  lies  in  the  range  7r/2^z,  <Tr.  An 
interesting  feature  of  this  system  is  the  growth  of  the  side 
wells  of  the  potential  as  the  nonlinearity  P  is  increased. 
This  feature  is  not  present  in  the  Josephson  junction, 
where  all  peaks  and  wells  of  the  potential  are  equal. 

In  Sec.  II  we  will  consider  the  time-dependent  analytic 
solution  of  (1.2)  on  the  separatrix.  This  is  accomplished 
via  a  spline  polynomial  approximation  of  the  nonlinear 
term  in  (1.2),  a  procedure  that  will  be  discussed  in  some 
detail.  This  method  has  also  been  carried  through  for  the 
off-separatrix  solution  and  is  discussed  briefly  in  this  sec¬ 
tion.  In  Sec.  Ill  we  apply  this  solution  to  the  calculation 
of  the  Melnikov  test  function  for  the  rf  SQUID  for  the 
center  as  well  as  the  side  wells  of  Fig.  1.  In  this  section 
we  also  numerically  search  for  the  homoclinic  points  and 
Cantor  set  structure  predicted  by  the  Melnikov  function, 
comparing  our  findings  with  the  theoretical  predictions. 
Finally,  in  Sec.  IV  we  display  evidence  for  a  strange  at¬ 
tracting  set  above  the  homoclinic  threshold.  This  numer¬ 
ical  result  is  discussed  briefly. 


II.  APPROXIMATE  ANALYTIC  SOLUTION 
OF  THE  rf  SQUID  EQUATION 

For  the  purpose  of  obtaining  a  separatrix  solution,  we 
now  turn  to  Eq.  (1.2)  and  formally  integrate  it,  obtaining 

'  P=vf+2mz,)-2mz)  ,  (2.1) 

where  (2,,u,-=z,-)  are  the  initial  values  used  to  determine 
the  integration  constant.  Equation  (2.1)  is  now  formally 


integrated  to  yield 

/(zi-/(.,).  I/;-/; 

dy 

[uf+2U{Zi)-2U{y)Y'^ 

=±t , 

(2.2) 

wnere  we  select  the  negative  sign  in  (2.2)  since  we  shall 
always  consider  z<z,  (this  choice  of  sign  yields  the 
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correct  monotonic  decreasing  behavior  on  the  separa- 
trix). 

In  order  to  evaluate  the  integral  in  (2.2),  we  employ  a 
polynomial  approximation  of  the  trignometric  function  in 
U (z).  Specifically,  we  set 

cosz^al—a2Z^+a^z*,  (2.3a) 


f ^(e)=vl+2U{Zi)—o)l 


2 

+4ircolpsm0  . 


(2.6d) 


One  readily  observes  that  for  any  0<z  <277,  the  integral 
in  (2.4)  takes  the  form 


where  02=0.4967  and  04=0.03705.  Such  approxima¬ 
tions  have  been  derived'^  by  expanding  cosz  in  a  contin¬ 
ued  fraction  series  that  can  be  truncated  at  any  order,  de¬ 
pending  on  the  level  of  accuracy  desired.  Differentiating 
(2.3a)  we  obtain 

sinzssojz— 03Z^  ,  (2.3b) 


where  O]  =0.9934  and  03=0.1482.  The  above  approxi¬ 
mation  (which  may  be  extended  to  include  higher-order 
terms  for  greater  accuracy)  has  been  estimated  to  yield  an 
error  of  0.13%  or  less  in  the  computation  of  the  trig¬ 
nometric  functions.  In  applying  these  approximations, 
one  is  under  the  restriction  z  <  v/2.  Hence,  for  z  >  ir/2, 
the  angle  z  must  be  expressed  in  terms  of  an  equivalent 
angle  in  the  first  quadrant  before  computing  the  trig¬ 
nometric  function.  We  shall  demonstrate  below  how  this 
spline  approximation  is  implemented. 

Let  us  assume,  as  an  example,  that  2v/2<z  <2it.  We 
break  the  first  integral  in  (2.2)  into  integrals  over  each  of 
the  quadrants  [a  similar  procedure  is  followed  for  the 
second  integral  in  (2.2)]: 


»/  V  f  ff/Z  rtr  riir/2 

1{Z)=  J  +  J-  +  J  + 

I  0  •'  ir/2  •'  <r 


X 


_ in. _ 

[vl+2mZi)-2my)Y^^  ■ 


(2.4) 


The  first  three  integrals  on  the  right-hand  side  of  (2.4) 
may  be  cast  in  the  form 


de 


1/2  ’ 


k  =  l,2,3 


(2.5) 


where  0<6< ir/2  so  that  the  expansions  (2.3)  are  valid. 
Finally,  the  last  integral  on  the  right-hand  side  of  (2.4) 
may  be  written  as  l4—I^(2ir—z),  where  takes  on  the 
form  of  (2.5)  (with  k  =4)  and  /^(z)  is  given,  in  general,  by 
an  integral  of  the  form  (2.5)  with  the  upper  limit  replaced 
by  z.  In  using  the  expansions  (2.3),  it  is  necessary  to 
transform  the  integrals  in  (2.4)  so  that  the  argument  of 
the  trignometric  function  is  restricted  to  the  first  qua¬ 
drant.  This  necessitates  the  breakup  of  the  integral  as 
above,  with  the  integrand  being  redefined  in  each  of  the 
quadrants  according  to 


fi{d)=v]+2U{Zi)—(ol 


2 

— 477(Uo/3sin0  , 


(2.6a) 


f2{e)=v^+2mZi)-(ol 


2 

— 477620)3 sin0 , 


(2.6b) 


/3(0)=v?+2U(z,)-<uI 


2 

-)-4776Joi3sin0  ,  (2.6c) 


f (z)~/ J  -r /2-i- ^3-^/4  —74(277 — z ), 


~<z<2v  (2.7a) 


—  I\-\-Ii-\-l3(z  — 77), 

2it 

77<Z  <  — 

-  -  2 

(2.7b) 

=  /,-1-/2-/2(77-z). 

77 

—  <Z  <77 

2  -  - 

(2.7c) 

=/,(z),  0<Z<^  . 

(2.7d) 

The  extension  of  the  above  procedure  for  z>2tt  is 
straightforward  and  will  not  be  discussed  here.  We  have 
thus  reduced  the  integral  (2.2)  to  a  sum  of  integrals,  each 
of  which  may  be  analytically  evaluated  using  the  approxi¬ 
mations  (2.3).  For  a  general  initial  and  final  value  of  z, 
the  integrals  I  are  evaluated  in  terms  of  the  Jacobian  el¬ 
liptic  functions.  On  the  separatrix,  however,  the  solution 
simplifies  considerably,  as  will  be  apparent  in  what  fol¬ 
lows.  We  might  mention  that  had  we  solved  the  dynamic 
equation  in  its  original  form  [Eq.  (l.D),  the  second  term 
in  (1.3)  would  have  contained  cosx,  in  which  case  we 
would  have  had  to  approximate  it  using  (2.3a).  The  func¬ 
tions  / k  appearing  in  (2.5)  and  (2.6)  would  then  be  quar- 
tics  rather  than  the  relatively  simple  cubics.  The  solution 
in  this  case,  while  still  analytically  tractable  (and  slightly 
more  accurate),  becomes  extremely  complicated;  in  par¬ 
ticular,  an  analytic  computation  of  the  Melnikov  function 
becomes  impossible.  This  is  the  reason  for  making  the 
phase  transformation  x—»-z  which  results  in  Eq.  (1.2). 

Before  proceeding  with  the  evaluation  of  the  solution 

(2.2) ,  we  briefly  consider  the  accuracy  of  the  potential 

(1.3)  in  light  of  the  approximation  (2.3b).  The  sinz  term 
in  (1.3)  is  replaced  by  sin(7r— z)  in  the  second  quadrant, 
— sin(z— 77)  in  the  third  quadrant,  and  — sin(27r— z)  in 
the  fourth  quadrant.  For  z  >  277,  the  procedure  is  repeat¬ 
ed.  The  expansion  (2.3b)  is  then  used  to  compute  U  (z). 
In  Table  I  we  list  the  values  of  the  turning  point  Zj  ob¬ 
tained  directly  from  (1.6)  and  through  use  of  the  approxi¬ 
mation  (2.3)  in  (1.6).  This  is  followed  by  a  computation 
of  the  values  of  the  potential  C/(Z|)  for  each  of  these 
values  of  Z|,  where  we  have  used  a  combination  of  (2.3) 
and  (1.3)  to  compute  the  approximate  values  of  f/(zj). 
Finally,  the  point  Z2  is  computed  via  the  condition 
t/(z,)=C/(z2),  using,  once  again,  the  appropriate  values 
of  C/(z,),  corresponding  to  the  direct  and  approximate 
[using  (2.3)]  calculations.  The  procedure  is  repeated  for 
different  P  values  with  (62o,Xo)=(l>0^  throughout.  This 
table  highlights  the  deviations  between  the  direct  and  ap¬ 
proximate  quantities,  with  the  error  introduced  in  the 
computation  of  the  turning  point  Zj  appearing  to  increase 
with  increasing  )3.  This  is  due  partly  to  the  error  inherent 
in  numerically  obtaining  the  turning  points  of  the  poten- 
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TABLE  I.  Parameters  relevant  to  the  separatrix  calculations  involving  the  potential  U(z)  and  its  approximation  using  the  cxpan- 
sion(2.3).  6)o=l,Xo=0. 


p 

Direct 

7, 

Approximate 

U(z^] 

Direct 

1 

Approximate 

Direct 

Zl 

Approximate 

0.8 

2.5223168 

2.498  5890 

11.294425  50 

11.29266308 

3.759  84842 

3.75639995 

0.9 

2.332681  1 

2.3105992 

11.710069  19 

11.71984893 

4.195  87066 

4.206957  30 

1.0 

2.2180747 

2.201  3499 

12.190044  30 

12.20925944 

4.479  62261 

4.48509690 

2.0 

1.8461779 

1.8693389 

17.93074521 

17.96423620 

5.60904139 

5.61104063 

2.5 

1.7861685 

1.8191063 

20.97966744 

20.991  m  30 

5.851  80765 

5.851  562  52 

tial  [through  solving  Eq.  (1.6)],  but  also  to  the  fact  that 
the  absolute  error  introduced  into  the  calculation 
through  the  use  of  the  approximations  (2.3)  increases 
with  increasing  nonlinearity. 

A.  Separatrix  solution 

We  now  turn  our  attention  to  the  formal  solution  (2.2) 
and  evaluate  it  explicitly  on  the  separatrix.  This  is  fol¬ 
lowed  by  an  example  using  a  representative  set  of  system 
parameters.  We  assume  that  the  particle  starts  at  the 
point  z,  (=72)  at  time  /=0  with  zero  initial  velocity.  It 
arrives  at  the  point  z,  at  t  =  00.  We  consider  the  case  of 
moderate  P  (  <  1.5)  for  which  the  points  z^  and  Z2  are 
one  quadrant  apart,  although,  as  will  be  indicated  later, 
the  extension  to  higher  P  values  is  readily  accomplished, 
albeit  somewhat  more  tediously.  Since  the  final  point  z  is 
always  contained  in  the  second  quadrant  (for  any  value  of 
P),  (2.2)  takes  the  form  (noting  that  rr  ^  Z;  <  3ir/2) 

Ii{2i—n)+l2{ir-z)=t  ,  (2.8) 

where  we  have  used  (2.7b)  and  (2.7c).  Consider  the  func¬ 
tion  fiiz)  appearing  in  the  integral  Iiiv—z).  This  func¬ 
tion  is,  generally,  a  cubic  in  z.  On  the  separatrix,  howev¬ 
er,  it  may  readily  be  seen  that  the  function  /2  has  only 
two  roots,  one  of  these  roots  being,  in  fact,  a  turning 
point.  In  other  words,  we  may  write,  on  the  separatrix, 

/2(z)=^2(2- !«!  I  )^(z-h  IojI  )  ,  (2.9) 

where  /42=4ir/?<u^3.  Here  one  sees  that  the  z  axis  is 
tangential  to  the  curve  /2(z)=0  at  the  point  z=  la||. 
In  writing  /2  in  the  form  (2.9),  the  location  of  the  roots 
on  the  7,  axis  is  determined  by  the  signs  appearing  in  the 
factors  on  the  right-hand  side.  We  also  have  (from  sim¬ 
ple  geometrical  considerations) 

Z/=7r-f  1  ay  I  =Z2  (2.10a) 

and 

z,=7r—  jci  I  .  (2.10b) 

The  integrals  in  (2.8)  are  now  readily  evaluated  to  yield 

z(t)=z,— atanh^^t  ,  (2.11) 

wherea=Z2— Zi  and 

Equation  (2.11)  is  the  solution  of  the  dynanuc  equation 
(1.2)  on  the  separatrix  (ZjZ,),  i.e.,  in  the  side  well  of  the 


potential.  One  sees  that  as  /— »±oo,  z—*zy,  and  the 
stable  and  unstable  manifolds  are  identical,  as  one  might 
expect.  Before  continuing  we  must  reiterate  the  fact  that 
the  special  form  of  the  solution  (2.1 1)  was  derived  for  the 
case  of  moderate  nonlinearity,  for  which  the  point  Z2  is 
contained  in  the  third  quadrant.  The  symmetry  between 
the  solutions  in  the  second  and  third  quadrants  [this  sym¬ 
metry  is  evident  through  a  glance  at  equations  (2.7b)  and 
(2.7c)]  allows  us  to  write  down  the  relatively  simple  ex¬ 
pression  (2.11).  The  foregoing  analysis  may  also  be  ap¬ 
plied  to  cases  (involving  larger  p  values)  for  which  the 
points  Zj  and  z^  are  several  quadrants  apart.  In  these 
cases,  however,  the  derivation  of  a  solution  analogous  to 
(2.11)  is  not  so  simple.  Specifically,  one  must  compute 
the  solution  z(t)  quadrant  by  quadrant  and  systematical¬ 
ly  piece  it  together  to  obtain  its  complete  behavior  for  all 
times.  The  solution  is  then  a  spline  function. 

Having  computed  the  separatrix  solution  (2.1 1)  it  is  in¬ 
structive  to  compare  it  with  the  solution  obtained  via  a 
direct  numerical  simulation  of  (1.2).  In  doing  so  we 
demonstrate  some  of  the  uncertainties  endemic  to  a  nu¬ 
merical  computation  of  the  separatrix  solution.  We  con¬ 
sider  as  a  specific  example  the  case  (^,(Uo,Xo)=(  1, 1,0). 
Writing  the  cubic  filz)  in  the  form 

2^^:)  =  A 27^ -i- B 27^ C 27 -{■  D 2  >  (2.12) 

we  obtain,  through  comparison  with  (2.6b), 

(^2.^2.^2,-02)^(1.862  34,-1,-3.058  65,  2.211  91)  , 

where  we  have  used  the  approximate  values  of  Zj,  Z2,  and 
(/(z  j )  from  Table  I  [since  the  approximation  (2.3)  is  con¬ 
tained  in  the  analytic  integration  of  (2.2)].  A  direct  nu¬ 
merical  computation  of  the  roots  of  the  cubic  (2.12)  yields 
three  roots,  two  of  which  lie  very  close  to  each  other. 
This  uncertainty  is  displayed  in  Fig.  2,  which  shows  a 
plot  of  the  function  /2(z).  The  point  z=ai  at  which  the 
curve  touches  the  z  axis  is  not  uniquely  determined, 
pointing  out  a  source  of  error  inherent  in  the  direct  nu¬ 
merical  computation  of  the  roots.  Such  a  numerical  com¬ 
putation  yields  the  roots  a  |  =0.94023-1-0.001 42/, 
02  =0.94023  —  0.00142/,  and  03=1.343  386.  It  is  ap¬ 
parent,  therefore,  that  even  though  the  approximate 
values  of  z,,  Z2,  and  I7(zi)  were  used  in  the  computation 
of  the  coefficients  of  the  cubic  (2.12),  numerically  solving 
the  cubic  yields  quantities  that  are  different  from  those 
that  were  input  into  the  coefficients  of  the  cubic  in  the 
first  place.  The  roots  Oj  and  03  coincide  if  one  ignores 
their  small  imaginary  parts  (the  existence  of  these  imagi¬ 
nary  parts  n.uy  be  directly  attributed  to  errors  inherent  in 


37 


HOMOCLINIC  CHAOS  IN  THE  rf  SUPERCONDUCTING  . . . 


3545 


FIG.  2.  Difference  U(z,)-U(z)  on  the  separatrix. 
(j3,wo,Xo)=(  1,1,0). 


the  numerical  computation  of  these  roots);  the  resulting 
root  represents  the  turning  point  at  which  the  curve  of 
Fig.  2  touches  the  z  axis.  In  Fig.  3  we  plot  the  solution 
zU)  on  the  separatrix  for  the  special  case  under  con¬ 
sideration  in  this  paragraph.  The  solid  curve  corresponds 
to  the  solution  obtained  via  a  direct  numerical  integra¬ 
tion  of  the  differential  equation  (1.2).  In  carrying  out  this 
integration,  the  direct  values  of  z^  and  Zj  used  are  from 
Table  I.  The  approximate  analytical  solution  (2.11)  is 
also  plotted  in  this  figure  (dotted  curve).  It  is  evident  that 
although  the  behavior  of  the  two  solution  curves  is  very 
close,  there  are  differences,  most  notably  at  long  times  (as 


FIG.  3.  Solution  zit)  [Eq.  (2.11)]  in  side  well  for  /?=!.  The 
dotted  curve  represents  Eq.  (2.11)  and  the  solid  curve  is  ob¬ 
tained  by  direct  (numerical)  integration  of  (1.2). 


one  would  expect  from  Table  I).  It  is  certainly  true  that 
in  addition  to  the  approximations  associated  with  the 
solution  (2.11),  tolerance-related  errors  may  be  intro¬ 
duced  via  the  numerical  simulation  of  Eq.  (2.1).  These 
latter  errors  are  more  difficult  to  estimate. 


B.  Off-separatrix  solution 


The  treatment  of  Sec.  II A  may  be  modified  and  ex¬ 
tended  to  cover  the  case  when  the  initial  value  z,-  does  not 
correspond  to  the  point  z^.  Once  again  we  utilize  the 
polynomial  approximation  (2.3)  in  evaluating  the  integral 
(2.2).  However,  the  initial  velocity  must  be  taken  to  be 
nonzero.  As  in  Sec.  II A  the  integral  in  (2.2)  is  broken  up 
into  the  sum  of  separate  integrals  over  the  four  qua¬ 
drants.  In  this  general  case,  however,  the  cubic  will 
have  three  distinct  roots  (two  of  which  may  be  complex 
conjugates  of  each  other).  Our  treatment  in  this  section 
is  brief  since  we  are  not  concerned  with  the  off-separatrix 
solution  in  the  rest  of  this  work.  The  procedure,  howev¬ 
er,  is  worth  outlining  since  it  may  be  of  practical  interest. 

Let  us  assume  that  the  initial  position  z,-  and  initial 
(nonzero)  velocity  u,-  aie  arbitrarily  chosen  at  (=0.  At 
any  later  time,  the  position  z(t)<Zj  may  be  written 
down,  using  (2.2),  (2.4),  and  (2.7),  as 

,  (2.13) 


where  I  takes  on  the  values  z,  n—z,  z—ir,  or  2it—z,  de¬ 
pending  on  the  quadrant  of  location  of  z  (the  cubic  ap¬ 
pearing  in  the  integrand  of  I  must  also  be  suitably  select¬ 
ed).  We  now  assume  that  the  cubic  f^iz)  has  three  dis¬ 
tinct  roots  «jti>ajfe2>a*3)  where  k  defines  the  quadrant 
of  location  of  the  point  z.  Then  one  may  write  (2.13)  in 
the  form 


dd 


0  [fkm 


1/2 


(2.14) 


where 


IAic2  being  the  coefficient  of  z^)  and  ^2  is  a  constant  to  be 
determined  by  the  initial  condition.  The  integral  /  (I)  is 
an  elliptic  integral  of  the  first  kind  and  one  may  cast 
(2.14)  in  the  form 


sn'(f,»+^2)  ’ 


(2.15) 


where  sn  is  the  elliptic  function  of  Jacobi  and  the  initial 
condition  enables  us  to  set  the  constant  ^2  via  the  condi¬ 
tion 


snf2= 


(2.16) 


In  practice,  it  is  often  more  convenient  to  evaluate  the 
Jacobian  elliptic  functions  in  terms  of  their  inverses,  i.e., 
to  evaluate  I(z)  directly  in  terms  of  the  elliptic  integrals 
of  the  first  kind.  Then  one  may  write  (2.15)  in  the 
equivalent  form, 
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!"»*,)]  .  (2.17) 

where  we  have  set 

<t>=  - - 

and  F  is  the  elliptic  integral  of  the  first  kind  having  am¬ 
plitude  m/^,  the  latter  quantity  being  expressed  in  terms 
of  the  roots  of  the  cubic  /^.  The  techniques  for  the  eval¬ 
uation  of  these  integrals  are  well  known'*  and  will  not  be 
repeated  here. 

III.  THE  MELNIKOV  FUNCTION 
FOR  THE  rf  SQUID 

We  now  suppose  that  the  unperturbed  system  dis¬ 
cussed  in  Sec.  II  is  perturbed  by  a  combination  of  dissipa¬ 
tive  and  periodic  forces,  i.e.,  F^  and  F^,  defined  in  Eqs. 
(1.4)  and  (1.5),  respectively.  The  purpose  of  this  section  is 
to  theoretically  investigate  the  condition  for  the  onset  of 
homoclinic  behavior  in  the  presence  of  the  above  pertur¬ 
bations. 

The  Melnikov  test  function  may  be  written  as,^“^’’ 

Al/olsf"  [F*(t)-i.F„(/)]i(/)dt  =  A*-t-A„(/o).  (3.1) 

This  remarkable  test  function,  valid  under  weak  perturba¬ 
tion,’  has  the  following  properties:  (a)  A(/o)=0  for  no 
perturbation  and  (b)  Altg)  changes  sign  (as  a  function  of 


tg),  indicating  a  crossing  of  the  stable  and  unstable  mani¬ 
folds,  i.e.,  the  presence  of  a  homoclinic  point.  As  dis¬ 
cussed  in  Sec.  I,  this  indicates  the  presence  of  a  Cantor 
set  structure  or  homoclinic  instability.’’*  The  simplicity 
of  this  procedure  arises  from  the  fact  that  the  quantity 
z(t)  in  (3.1)  |s  the  unperturbed  separatrix  velocity  in  the 
region  (ZjZj),  given  approximately  by  the  time  derivative 
of  Eq.  (2.11).  This  arises  because  the  separatrix  orbits 
are,  in  a  sense,  the  most  sensitive  to  any  perturbation. 
The  Melnikov  function  can  only  be  applied  to  separatrix 
orbits.  No  analysis  has  been  carried  out  for  nonsepara- 
trix  orbits  and  thus  it  is  applicable  only  to  unperturbed 
nonlinear  oscillators  with  a  hyperbolic  fixed  point  and  a 
separatrix.  The  rf  SQUID  is  a  good  example.  In  apply¬ 
ing  the  Melnikov  function  to  the  rf  SQUID,  considerable 
care  must  be  taken  to  ensure  that  the  magnitudes  of  the 
parameters  P,k,  A,  etc.  or  combinations  of  these  parame¬ 
ters  fall  within  the  realm  of  validity  of  perturbation 
theory. 

Let  us  now  evaluate  the  integrals  in  (3.1),  considering 
first  the  region  (zjZj)  of  the  potential  (1.3)  (see  Fig.  1). 
Using  the  solution  (2.11)  (we  are,  once  again,  confining 
ourselves  to  the  moderate  )3  case),  we  readily  obtain  the 
first  integral  as 

A*  =  |fa’^k,  (3.2) 

where  the  quantities  a  and  ^  have  been  defined  in  connec¬ 
tion  with  Eq.  (2.11).  In  order  to  evaluate  the  second  in¬ 
tegral  in  (3.1)  we  express  it  as  the  sum  of  two  integrals: 


AJto)  =  —2a^Af°°  tanh^t  secb^^tisinot  cosatQ—costot  sinci)to)dt  .  (3.3) 

—  00 


Noting  that  z(t)  is  an  odd  function  in  the  side  well  of  the 
potential,  the  above  integral  becomes 

—  4af  y4(cos<u/o)  fj‘  tanhft  sech’^/  sinwt  dt  . 

(3.4) 

This  integral  is  readily  evaluated  to  give 
4it(o^A  COStUtn 

AJ/o)=-  '  .  .  r  ,,,  ,1/21  »  (3.5) 

A2Sinh[Trco/(A2a]'] 

where  the  quantity  A^  has  been  defined  in  connection 
with  Eq.  (2.9).  It  is  evident  that  the  function  A„  has  its 
extrema  for  cosa)to=±l  [all  other  quantities  in  (3.5)  be¬ 
ing  fixed].  In  the  presence  of  finite  damping  one  obtains 
homoclinic  behavior  above  a  critical  threshold,  which 
may  be  found  by  setting  Ajt=  [  A^lrgl/costuto  I  ■  This 
condition  leads  one  to  the  threshold  condition  for  the  on¬ 
set  of  homoclinic  behavior, 

-r  =  — - r^sinh[7r<a/(/l2«)‘''^]  •  (3-6) 

K  lOTT  (o‘- 

A  similar  result  was  fir^t  obtained  by  Holmes'*  for  the 
anti-Duffing  oscillator.  Damping  may  suppress  chaos;  in 
this  case,  the  forward  manifold  spirals  from  the  hyperbol¬ 
ic  fixed  point  into  the  elliptic  fixed  point,  while  the  back¬ 


ward  manifold  spirals  outward;  the  manifolds  do  not 
cross.  As  mentioned  in  Sec.  I,  an  analysis  similar  to  that 
just  carried  out  has  been  done  for  the  simple  Josephson 
junction.’"'®  The  analysis  in  that  case  is  considerably 
simplified  by  the  fact  that  the  term  linear  in  x  in  Eq.  (1.1) 
is  absent  from  the  Josephson-junction  dynamics,  and  the 
unperturbed  solution  may  be  obtained  analytically 
without  the  approximations  of  Sec.  II.  Unlike  the  rf 
SQUID,  the  Josephson  junction  is  described  by  a  poten¬ 
tial  in  which  all  the  wells  have  the  same  size;  in  fact,  the 
potential  in  this  case  is  simply  of  the  form  t/(z)~sinz. 
For  the  system  under  consideration  in  this  work,  one  re¬ 
covers,  qualitatively,  the  features  of  the  Josephson  junc¬ 
tion  in  the  limit  /?->  oo  (in  this  limit,  all  the  wells  in  the 
potential  of  Fig.  1  approach  the  same  size).  In  what  fol¬ 
lows  we  shall  compare  the  results  of  this  section  with  ex¬ 
isting  results  for  the  Josephson  junction. 

In  Fig.  4  we  plot  the  quantity  i  A^tg^/costa/o  |  as  a 
function  of  the  driving  frequency  to  for  li=  1  in  the  side 
well  (Z|  <z<Z2)  of  Fig.  1.  The  natural  frequency  (Oq  is 
set  equal  to  unity  in  this  plot.  The  constant  quantity  A^ 
of  Eq.  (3.2)  is  also  shown  (the  straight  line).  The  solid 
curve  is  obtained  by  numerical  evaluation  of  the  integral 
in  (2.2)  and  a  subsequent  numerical  computation  of  the 
Melnikov  integral  (3.4).  Also  shown  are  data  points  cor¬ 
responding  to  the  theoretical  result  obtained  from  (3.5j. 
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FIG.  4.  Melnikov  function  amplitude  in  the  side  well  (Z1Z2). 
The  solid  curve  represents  numerical  calculation  and  data 
points  are  derived  from  Eq.  (3.5).  The  dashed  line  represents 
Eq.  (3.2).  (P,k,A)={l,lA). 


Note  that  the  difference  between  a  numerical  calculation 
of  Aj^  and  the  theoretical  result  (3.2)  is  much  smaller  than 
the  scale  of  this  figure.  For  the  moderate  value  of  p  used 
in  this  figure,  the  agreement  between  theoretical  and  nu¬ 
merical  results  is  quite  good.  We  reiterate  that  the  nu¬ 
merical  computations  are  not  totally  accurate;  however, 
it  is  difficult  to  obtain  realistic  estimates  of  the  error  in¬ 
troduced  into  the  numerical  integrations  of  Eqs.  (2.2)  and 
(3.4).  The  intersection  of  the  straight  line  with  the  curve 
represents  the  threshold  for  the  onset  of  homoclinic  be¬ 
havior  since  the  Melnikov  function  changes  sign  as  one 
crosses  the  line.  The  graph  of  j  A^(to)/cos<uro  |  is 
peaked  at  m=<Uc,  given  [for  the  theoretical  result  (3.5]  by 


This  peak  represents  the  minimum  value  oi  ( A /k)  neces¬ 
sary  for  the  onset  of  homoclinic  behavior.  For  P—  1  one 
obtains  the  critical  values  (a),/4A)=(  1.257,1.79)  from 
the  theoretical  calculations.  These  numbers  compare 
quite  favorably  with  the  values  (1.26,1.78)  obtained  via 
the  numerical  calculations. 

Figure  5  shows  the  quantity  |  A^(to)/sin«to  |  for  the 
center  well  of  the  potential  (1.3)  (i.e.,  — tt— z,<2<z,) 
corresponding  to  the  separatrix  (Z1Z3)  of  Fig.  1.  For  the 
rf  SQUID  we  expect  the  results  in  the  center  well  to  differ 
from  those  obtained  in  the  side  well,  both  because  of  the 
difference  in  well  dimensions  and  the  difference  in  the 
parity  of  the  velocity  between  the  wells  (both  these 
features  are  absent  in  the  Josephson  junction).  This 
figure  has  been  obtained  numerically.  If  one  considers 
the  velocity  profiles  in  the  two  wells  (Fig.  6),  it  is  seen 
that  the  velocity  in  the  side  well  is  odd.  Hence,  in  (3.3), 
only  the  first  term  contributes  to  the  integral  which  is 
peaked  as  some  nonzero  frequency  this  corresponds 
to  the  case  worked  out  analytically  in  this  section.  In  the 
center  well,  however,  the  velocity  is  even  so  that  the  Mel¬ 
nikov  integral  in  this  case  is  proportional  to  sincatQ  in 
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FIG.  5.  Melnikov  function  amplitude  in  the  center  well 
(zjZj)  for  the  same  parameter  values  as  Fig.  5.  Both  curves  are 
obtained  numerically. 


agreement  with  results  derived  for  the  Josephson  junc¬ 
tion.^"**  In  Fig.  7  wc  show  the  critical  (i.e.,  minimum) 
value  ( i4  A )(.  as  a  function  of  the  nonlinearity  parameter 
P,  where  the  subscript  c  implies  that  we  have  evaluated 
( A  /k)  at  the  critical  frequency  corresponding  to  the 
side  well.  Results  for  the  center  as  well  as  the  side  well  of 
the  potential,  are  plotted.  As  expected,  when  a  greater 
disparity  exists  in  the  relative  dimensions  of  the  wells  (for 
moderate  p),  a  higher  minimum  value  of  {A/k)  is  re¬ 
quired  to  trigger  homoclinic  behavior  in  the  center  well 
at  the  critical  driving  frequency  0)^  corresponding  to  the 
side  well.  For  large  P  values,  the  values  of  ( i4  //c)^  in  the 
two  wells  converge. 


FIG.  6.  Velocity  profile  i(i)  in  the  center  well  (solid  curve) 
and  in  the  side  well  (dotted  curve).  P=  1. 
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FIG.  7.  Critical  values  of  ( .4  //c)  as  a  function  of  /3  in  center 
well  (solid  curve)  and  side  well  (dotted  curve).  Both  curves  are 
derived  numerically  at  the  critical  frequency  corresponding 
to  the  peak  of  the  Melnikov  function  amplitude  in  the  side  well. 


IV.  HOMOCLINIC  THRESHOLD 
AND  “STRANGE”  ATTRACTOR  (NUMERICAL) 

In  this  section  we  will  show  Poincare  return  maps  of 
the  numerical  solutions  to  the  driven  damped  rf  SQUID, 
Eqs.  (1.2)-(1.5).  The  objective  will  be  to  first  estimate  the 
onset  of  homoclinic  crossing  and  compare  the  results 
with  those  of  Sec.  III.  This  was  first  done  by  Holmes^  for 
a  cubic  map,  which  preserved  some  of  the  properties  of 
the  anti-Dufiing  oscillator.  No  such  results  have  previ¬ 
ously  been  obtained  for  the  Josephson  junction  or  the  rf 
SQUID.  In  the  final  part  of  this  section  we  will  give  evi¬ 
dence  for  the  existence  of  a  strange  attracting  set  in  the 
solutions,  which  develops  in  a  parameter  range  beyond 
the  onset  of  the  homoclinic  instability.’  Such  attractors 
have  been  numerically  studied  in  detail  for  the  Josephson 
junction  (see  Ref.  8  for  a  review). 

Comments  are  made  by  Holmes'*  concerning  numerical 
techniques,  but  it  is  worthwhile  to  also  comment  here  on 
the  procedure  and  the  difficulties.  In  order  to  facilitate 
the  numerical  integration  of  the  differential  equation  (1.1) 
we  have  introduced  the  scaled  time  variable  T=aQt.  It  is 
apparent  that  for  the  special  case  cjq  considered 
throughout  this  work,  this  scaling  does  not  change  the 
original  equation  at  all.  However,  in  all  practical  appli¬ 
cations  of  rf  SQUID’S,  one  usually  has  lO'®,  which 
leads  to  enormous  problems  when  one  attempts  to  in¬ 
tegrate  the  equation  of  motion  (1.1),  unless  such  a  scaling 
is  utilized.  The  stable  and  unstable  manifolds  discussed 
in  this  section  are  computed  by  mapping  a  large  number 
of  points  on  the  stable  (unstable)  manifold  near  the  saddle 
point,  one  or  more  Poincare  periods  backward  (forward). 
In  the  presence  of  damping  and  the  periodic  external  per¬ 


turbation  one  is  faced  with  the  problem  of  computing  the 
saddle  point.  Although  numerous  methods  exist  for  this 
purpose  (e.g.,  averaging’'*®  and  harmonic  balance*®),  the 
answer  obtained  for  the  form  of  the  nonlinearity  in  this 
problem  is  not  accurate  enough.  Nevertheless,  by  es¬ 
timating  the  location  of  the  saddle  point  and  mapping  a 
small  area  around  this  point  forward  and  backwards  in 
time  [using  Eq.  (1.1)],  one  obtains  a  reas^'nably  good  pic¬ 
ture  of  these  manifolds  and  their  behavior  relative  to 
each  other.  In  order  to  obtain  clear  and  distinct  pictures 
of  the  manifolds,  one  must  make  the  mapped  area  (about 
the  saddle  point)  arbitrarily  small.  In  the  system  at  hand, 
however,  one  cannot  do  this  since  the  saddle  point  has 
not  been  accurately  located.  This  uncertainty  is  most 
likely  the  reason  for  the  smearing  of  the  manifolds  in 
Figs.  8-13.  It  is  also  possible  that  some  numerical  uncer¬ 
tainties  arise  because  it  is  difficult  for  the  integrating  rou¬ 
tine  to  exactly  follow  the  true  solution  of  the  differential 
equation  so  close  to  an  unstable  fixed  point. 

Let  us  now  consider  some  results.  We  will  take 
()3,ft),A:)=(2,2.25, 1)  (the  value  (0=2.25  is  chosen  be¬ 
cause  it  is  very  close  to  critical  frequency  (o^  for  this 
value  of  the  nonlinearity)  and  vary  the  driving  amplitude. 
We  set  q  =  A  /2it  and  note  that  throughout  this  section 
we  work  in  the  original  system  of  variables  (x,x,t)  of  Eq. 
(1.1).  For  q=0.5  shown  in  Fig.  8,  the  unstable  hyperbol¬ 
ic  fixed  point  is  close  to  the  'unperturbed  hyperbolic  point 
.X  =0.5438  (this  point  corresponds  to  the  point  Zy  in 
Table  I).  A  branch  of  ine  unstable  manifold  spirals  into 
the  right-hand  elliptic  fixed  point  due  to  the  dispersion 
{k  =  \).  A  branch  of  the  stable  manifold  emanates  from 
the  region  of  the  center  well  and  approaches  the  hyper¬ 
bolic  point  “around”  the  unstable  solution  in  the  Poin¬ 
care  phase  plot;  the  stable  and  unstable  manifolds  do  not 
touch.  This  pattern  is  characteristic  of  an  overdamped 
system.  For  q= 0.72  shown  in  Fig.  9,  the  unstable  mani¬ 
fold  approaches  the  stable  manifold  by  developing  a 
sharp  cusp.  At  q  =0.73  (Fig.  10)  the  manifolds  appear  to 
touch,  and  in  Fig.  11  they  have  already  crossed,  for 


FIG.  8.  Stable  and  unstable  manifolds  in  parameter  legime 
where  no  homoclimc  behavior  is  expected.  i(},(o,k,q  =  A/ 
25r)=(2,2.25, 1,0.5). 
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x(t) 

FIG.  9.  Same  as  Fig.  8  but  with  9  =0.72.  The  stable  mani¬ 
fold  has  developed  a  cusp  and  is  approaching  the  unstable  mani¬ 
fold  in  the  side  well. 


q  =0.)5,  in  the  side  well.  For  these  parameter  values,  the 
Melnikov  function  of  Sec.  Ill  predicts  a  homoclinic 
crossing  for  9  =0.74,  in  good  agreement  (considering  the 
numerical  difficulties)  with  our  numerical  results.  Fur¬ 
ther  increasing  q  causes  the  other  branch  of  the  stable 
manifold  to  loop  back  and  cross  the  other  branch  of  the 
unstable  manifold;  the  latter  is  attracted  to  the  stable  (el¬ 
liptic)  fixed  point  in  the  center  well.  Figure  12  shows  the 
near  crossing  for  9  =  1.21,  and  in  Fig.  13  the  homoclinic 
intersection  has  already  occurred  (for  9  =  1.25).  The 
Melnikov  function  predicts  the  first  crossing  in  the  center 
well  for  these  parameters  at  9  =  1.2. 

Now  let  us  turn  to  the  evidence  for  the  appearance  of  a 
global  steady  chaotic  attracting  set  in  the  SQUID — a 
strange  attractor.  Chaos  in  such  a  set  may  be  properly 
characterized  by  its  Lyapunov  exponent. In  Fig.  14 
we  plot  the  maximal  Lyapunov  exponent  as  a  func¬ 
tion  of  the  periodic  driving  amplitude  q  with 
{P,(OQ,k)=(2,l,l).  This  quantity  was  computed  using 


FIG.  10.  Same  as  Fig.  9  but  with  q  =  0.1Z.  Critical  case,  the 
manifolds  just  touch  in  the  side  well. 


x(t) 

FIG.  11.  Same  as  Fig.  10  but  with  9=0.75.  Supercritical 
case;  a  homochnic  crossing  has  taken  place  in  the  side  well. 


the  algorithm  of  Wolf  ct  a/.^'  It  is  evident  that  the  system 
displays  chaotic  behavior  (characterized  by  a  positive  A.„) 
at  a  value  of  q  somewhat  above  the  homoclinic  threshold 
value  9  =  1.21.  Further,  one  observes  bands  of  periodic 
behavior  (characterized  by  a  negative  A.,„)  at  higher  q 
values.  Such  intermittent  behavior  is  now  known  to 
occur  in  many  nonlinear  chaotic  systems  and,  in  particu¬ 
lar,  in  the  driven  Duffing  oscillator.'* 

We  should  comment  that  there  is  no  theoretical  con¬ 
nection  between  to  the  appearance  of  these  attracting  sets 
and  what  we  have  termed  homoclinic  instability.  The 
Melnikov  function  does  not  characterize  the  appearance 
of  a  strange  attractor,  as  is  well  known.*  There  have  been 
efforts  to  correlate  the  Melnikov  test  empirically  with  the 


x(t) 

FIG.  12.  Same  as  Fig.  11  but  with  9  =  1.21.  Near-critical 
case  for  homoclinic  crossing  in  the  center  well;  the  manifolds 
are  about  to  touch  in  the  center  well  (upper  left  of  figure).  For 
this  value  of  9,  homoclinic  crossing  has  already  occurred  in  the 
side  well.  The  blackened  area  represents  the  smearing  of  the 
manifold  due  to  the  numerical  uncertainties  referred  to  in  the 
text. 
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FIG.  13.  Same  as  Fig.  12  but  with  9  =  1.25.  Supercritical 
case;  homoclinic  crossing  has  just  occurred  in  center  well. 


appearance  of  “chaos”  in  the  Josephson  junction.*"" 
While  not  uninteresting,  they  do  not  evidence  a  physical 
prelude  or  early  scenario  to  the  appearance  of  a  strange 
attractor.  At  best,  now,  one  would  expect  that  the  Melni¬ 
kov  function  is  a  “rule  of  thumb”  insofar  as  the  appear¬ 
ance  of  a  strange  attractor  is  concerned.  It  provides  one 
with  a  lower  bound  for  the  chaos  threshold  in  a  given 
nonlinear  system  and  its  vanishing  should  be  considered  a 
necessary  condition  for  the  appearance  of  chaos  in  the 
dynamics  under  consideration.  Moon  and  Li^^  have  con¬ 
structed  the  fractal  basin  boundary  for  the  driven  anti- 
Duffing  oscillator  and  have  observed  that  the  fractal 
structure  appears  to  be  correlated  with  the  appearance  of 
homoclinic  orbits  in  this  system.  They  observe  that 
above  the  homoclinic  threshold  (determined  by  Holmes'* 
using  the  Melnikov  integral),  the  fractal  basin  boundary 
becomes  quite  complicated,  whereas  it  is  smooth  and 
nonfractal  below  this  threshold.  They  conclude  that  the 
Melnikov  criterion  is  a  necessary  condition  for  the  ap¬ 
pearance  of  the  complicated  fractal  boundary.  Similar 
boundaries  have  been  constructed  for  the  driven  damped 


FIG.  14.  Maximal  Lyapunov  exponent  as  a  function  of 
the  driving  force  amplitude  q.  (/?, to, /c  )  =  (2,2.25, 1 ). 


FIG.  15.  Attractor  corresponding  to  9  =  1.43,  jfc  =  1.0.  This 
case  represents  the  possible  onset  of  a  chaotic  attractor. 


pendulum  (this  is  the  mechanical  analog  of  the  Josephson 
junction)  by  Gwynn  and  Westervelt.^* 

We  now  examine  the  parameter  range  for  which  a 
chaotic  steady  attractor  seemingly  occurs.  At  any  given 
time,  the  state  of  the  system  is  completely  specified 
through  a  determination  of  its  position  and  velocity 
(x,x).  It  is  well  known  that  in  the  chaotic  regime,  a 
Poincare  plot  of  the  system  evolution  displays  a  strange 
attractor,  i.e.,  there  exist  steady  areas  in  state  space  to 
which  the  states  of  the  system  are  preferentially  attracted 
in  a  seemingly  random  manner — the  successive  values  of 
(x,x )  jump  from  one  region  of  state  space  to  another  in  a 
random  manner  producing  a  topologically  complex  Poin¬ 
care  section. 

In  Fig.  15  we  show  the  Poincare  plot  of  the  system  for 
9  =  1 . 43,  k  =  1 . 0;  Fig.  16  shows  the  effects  of  reducing  the 
damping  to  k  =0.3  (the  attractor  displays  a  far  more  in¬ 
tricate  structure  at  this  lower  k  value).  The  value 
9  =  1.43  is  seen  (from  Fig.  14)  to  be  quite  close  to  the 
threshold  for  the  onset  of  chaos  in  this  system.  The 
(common  logarithms  of  the)  power  spectral  densities  cor- 
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FIG.  16.  Chaotic  attractor  corresponding  to  0=1.43, 
A- =0.3. 
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FIG.  17.  Spectral  densities  corresponding  to  the  attractors  of 
Figs.  15  (solid  curve)  and  16  (dashed  curve).  Only  the  harmon¬ 
ics  of  the  driving  frequency  arc  present,  with  the  prominent 
peak  on  the  left  occurring  at  the  fundamental  driving  frequency 
w/2ir=0.36. 


responding  to  these  cases  are  plotted  in  Fig.  17,  with  the 
power  spectral  density  defined  by 


(r)exp(<Ar  )dt 


(4.1) 


where  T-*  oo  and  the  angular  brackets  denote  the  time 
average  (the  spectra  were  computed  over  KXX)  rf  cycles 
with  a  frequency  resolution  of  approximately  0.021  Hz). 
We  observe  that  the  power  spectra  contain  harmonics  of 
the  driving  frequency  {<u/2ir=0.36)  only,  superimposed 
on  a  broadband  noise  background;  this  is  a  characteristic 
of  chaotic  behavior.  Finally,  in  Fig.  18,  we  show  a  case 
in  which  the  system  is  driven  enough  strongly  so  that  it 
makes  excursions  to  numerous  side  wells  in  the  poten¬ 
tial.  The  parameters  used  in  this  case  are  {p,a,k,q) 
=(15,2.25,1.5,21).  The  attractor  displays  a  quasi- 
periodic  overall  structure. 

In  concluding,  it  is  worthwhile  to  speculate  on  the 
significance  of  the  frequency  defined  in  (3.7).  From 
(3.6)  one  observes  that  for  a  given  damping  k,  the  period¬ 
ic  force  amplitude  A  necessary  to  trigger  homoclinic  be¬ 
havior  is  a  minimum  for  this  value  of  the  periodic  driving 
frequency.  Huberman  et  al}*  have  examined  the  appear¬ 
ance  of  chaos  in  the  sinusoidally  driven  Josephson  junc¬ 
tion  as  a  function  of  the  driving  frequency  <o  of  the  exter¬ 
nal  signal.  They  observe  chaotic  behavior  over  a  band  of 
driving  frequencies  centered  about  <o^coj/2,  where  coj  is 


X(t) 


FIG.  18.  Chaotic  attractor  for  (/3,<u,/c,9)=(15,2.25,l. 5, 21). 
The  system  makes  excursions  to  numerous  side  wells. 


the  plasma  frequency  of  the  junction.  Ritala  and 
Salomaa'^  have  observed  the  appearance  of  subharmonic 
and  chaotic  solutions  in  the  rf  SQUID  in  a  band  of  fre¬ 
quencies  centered  about  <uw(u,/2.  Here,  tu,  is  the  reso¬ 
nance  frequency  that  governs  low-amplitude  oscillations 
in  the  rf  SQUID  and  is  easily  shown  to  be  given  by 
ti)J=6)o(  I  +2'np).  The  quantity  u),  may  be  considered  the 
analog  of  the  plasma  frequency  for  the  rf  SQUID.  It  is 
interesting  to  note  that  as  p  is  varied,  the  ratio  (0^/0,  ap¬ 
pears  to  lie  in  the  interval  [0.4,0.9]  for  moderate  0 
(0.74<)3^  1(X)).  Hence,  one  might  expect  that  setting 
the  driving  frequency  0  approximately  equal  to  the  op¬ 
timum  value  0f  predicted  by  the  Melnikov  function 
might  provide  a  connection  between  our  results  and  ear¬ 
lier  work.*’’**  Of  course  it  must  be  remembered  that  sim¬ 
ply  selecting  an  appropriate  driving  frequency  according 
to  the  above  “prescription”  is  not,  by  itself,  sufficient  to 
induce  homoclinic  behavior  or  chaos  in  the  system;  the 
values  of  the  damping  coefficient  k  and  the  driving  force 
amplitude  A  must  also  be  appropriately  selected. 
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